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We analyze the Physics of cold atoms in honeycomb optical lattices with on-site repulsion and 
spin-orbit couplings that break time reversal symmetry. Such systems, at half filling and large 
on-site repulsion, have been proposed as a possible realization of the Kitaev model. The spin-orbit 
couplings break the spin degeneracy and, if strong-enough, lead to four non-overlapping bands in the 
non-interacting limit. These bands carry non-zero Chern number and therefore the non-interacting 
system has non-zero angular momentum and chiral edge states at 1/4 and 3/4 filling. We have 
investigated the effect of interactions using the variational cluster perturbation theory and conclude 
that the chiral edge states exist in finite range of interaction and hopping parameter space. 
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INTRODUCTION 

Systems with topological quantum order have been of 
much interest recently. These systems are characterized 
by quasi-particles with fractional charge and statistics. 
Systems with non-Abelian statistics are especially inter- 
esting since they can be used for fault-tolerant quantum 
computation. The Kitaev model on the honeycomb lat- 
tice is a quantum spin-^ model that has such topologi- 
cal order pQ. There have been recent proposals [2H4] to 
realize this model in an optical lattice, using two- state 
bosonic or fermionic atoms with spin-dependent hopping 
on a honeycomb lattice with on-site repulsion. The Ki- 
taev model is then the effective Hamiltonian at half-filling 
and large repulsion. 

In this paper, we concentrate on the fermionic model 
defined on a honeycomb lattice. Let b a (a = 1,2,3) de- 
note the vectors linking a site on the A sublattice with its 
three neighbors on the B sublattice, and let (ij) a denote 
the nearest-neighbor pairs along these three directions. 
The Hamiltonian of the model is then 

H=Y, {4 pa Cj + H -c} + U X] n zt n zi (1) 

where P a = \(t + t f a a ) and a a are the Pauli matrices 
(we have suppressed the spin indices). The number of 
fermions at site i is rii a = c\ a Ci a . At t' = 0, the model re- 
duces to the simple spin- invariant, nearest-neighbor (NN) 
Hubbard model, as relevant to graphene, and is time- 
reversal invariant. At t' = t, the one-body part of the 
Hamiltonian is a combination of the projection operators 
7^(1 + <r a ). Thus, only electrons that are spin polarized in 
the a direction can hop along the a bonds. At this value 
of t', at half filling and large U, the effective low energy 
spin model is the Kitaev honeycomb model [2H4]. 

The term proportional to t' is a spin-orbit coupling. It 
is qualitatively different from the usual spin-orbit term 



as it breaks time reversal symmetry. The physical origin 
of this lies in the asymmetry of the laser beams coupling 
the two spin components to an excited level (which makes 
the potential barrier spin dependent). It is interesting to 
note that a time reversal invariant hopping term of the 
form H = Yl(ij) a (h c l + it'<j a ) Cj + H.c^ does not lead 
to the Kitaev model at half filling in the large- U limit. 

Topological effects in electronic bands were highlighted 
in the seminal paper of Thouless et al. [5 where the 
quantized Hall conductivity was expressed as the sum 
of the Chern numbers of the occupied bands (when the 
Fermi level lies in a gap). The Chern number and hence 
the quantized Hall conductivity was later identified with 
the number of chiral edge channels at the Fermi level by 
Hatsugai [6 . A non-zero Chern number can only occur 
in a system which is not time reversal invariant. How- 
ever, this does not necessitate an external magnetic field: 
Haldane constructed a tight-binding model on a honey- 
comb lattice with next-NN hopping terms that are not 
time reversal invariant [7 j. In that model the two bands 
carry Chern numbers equal to ±1 and an anomalous Hall 
effect occurs when the bands are partially filled [8 . The 
Hall conductivity is quantized when the Fermi level lies 
in the gap. More recently, the topology of time reversal 
invariant electronic bands with spin-orbit couplings has 
been extensively studied in the physical context of the 
spin Hall effect, leading to the discovery of topological 
insulators [9l HQ] , The effect of time reversal breaking 
spin-orbit terms, to our knowledge, has not been studied 
so far. 

In the rest of the paper, we concentrate on Model ([!]) 
at \ (or |) filling and show that it is in an anomalous 
quantized Hall state in a region of the (t'-U) plane. At 
U = 0, the four bands carry Chern numbers equal to 
±1. At large enough t', they do not overlap, leading to 
quantized anomalous Hall states at fillings \ and |. We 
analyze the corresponding chiral edge state structure and 
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propose a way to detect them. We then study the model 
at U > using Cluster Perturbation Theory [TQ (CPT) 
and the Variational Cluster Approximation [T2] . 

The model is easily diagonalized in the U = limit. 
See Supplemental Material for details. The single- 
particle Hamiltonian in momentum space can be written 
as 

ft(k) = + ® S(k) + a~ <g> E*(k)) (2) 

k 

where E(k) = P3 + Pie* fel + P 2 e~* /C2 , fci } 2 are the com- 
ponents of the lattice basis vectors ei 5 2, and a ± are the 
raising and lowering Pauli matrices acting on sublattice 
indices (A and B). The single-particle spectrum is com- 
pletely determined in terms of the spectrum of the 2x2 
positive semi-definite matrix, EE''' = Cq + \ {t') 2 + |-B • a 
where e (k) = ||(l + e i/ci + e~ i/C2 )| and 

B = (1 — t' sin ki + cos &2 + cos ^3) bi 
+ (1 + cos k\ — t' sin ki + cos £3) b 2 
+ (1 + cos k\ + cos &2 — sin £3) b3 (3) 

The energy spectrum of the four bands is e(k) pp / = 
pf (eg(k) + §(i') 2 +p|*'|B|) = , where p.j/ = ±1. 




FIG. 1. (Color online) The two positive energy bands at 
t' = 0.1 (top) and t' = 0.9 (bottom). The negative energy 
bands are symmetrically located underneath the zero-energy 
plane and do not touch the upper bands. The two bands are 
overlapping at t' — 0.1 and well separated at t' — 0.9 

At t' = 0, we recover the graphene spectrum: The are 
two spin degenerate bands touching each other at two 
distinct Dirac points at the zone boundary. At any non- 
zero value of t', the spin degeneracy is broken. Two of 

the bands, e + _and e , continue to touch each other at 

two Dirac points whereas the other two develop a gap. 

The e ++ , and e |_, e__ bands overlap until a critical 

value t' c of t' given by - = 0.717. For t' > t f c , a 
finite gap appears. This is illustrated in Fig. [l] 



At t' = the middle two bands, e + _ and e , touch 

each other at k = 0. For t' > \/3, the point of contact 
splits into six more Dirac points, for a total of eight Dirac 
points. 

The Chern number for each band is ±z/, where v is 
the winding number of the unit vector field B. It is not 
difficult to prove that |B| ^ for all points in the Bril- 
louin zone, for any non-zero value of t' and hence v is 
well defined. We have numerically computed v to be 1. 
The highest and lowest energy bands with energies, =be + , 
have Chern numbers +1 whereas the middle two bands 
with energies, ±e_, have Chern numbers —1. 

When t' > 0.717 and the bands are non-overlapping, 
the half-filled state has total Chern number zero but the 
quarter and three-quarter filled states have Chern num- 
bers equal to ±1 respectively. These states are thus anal- 
ogous to quantum Hall states with quantized Hall con- 
ductivity an = ±1. At t' < 0.717, the bands overlap and 
at quarter filling the Fermi level passes through the e_± 
bands. The system is then analogous to an anomalous 
Hall conductor with the "Hall conductivity" equal to the 
Pancharatnam-Berry curvature integrated over the occu- 
pied states |B]. 




FIG. 2. The spectrum for an open tube of circumference 
L = 60 with zig-zag edges. Note the edge states that cross 
the gaps between the bands. 

In the regime t' > 0.717, we expect chiral edge states 
in the gap between the e + ± bands and also in the gap 
between the e_± bands. Our numerical calculations con- 
firm this. Figure [2] shows the spectrum of the system 
with a cylindrical geometry and zig-zag edges. There are 
non-dispersive, zero energy edge states, as in graphene. 
There are also chiral edge states between the top two and 
the bottom two bands. 

In a system of neutral atoms, the Hall conductivity is 
not easily measurable. However, the non- trivial topology 
also manifests itself in the orbital angular momentum, 
which is easier to measure then. The orbital magnetiza- 
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tion of Bloch electrons of the band pp' is given by [T3HT7] 

d 2 k 

/(k)<„ ( 27r ) 2 

x {8*$™- \(H k +pe p ,(k) - 2p)\d k &*'), (4) 

where \i is the chemical potential (Fermi energy) and 
\§pp>) are the single particle eigenvectors. In the metallic 
case, Eq. Q provides a /i-dependent magnetization, as 
it should. In the insulating case, when fi is varied in 
the gap, M changes linearly only if the Chern invariant 
is non-zero, and remains constant otherwise. Eq. Q is 
related - though not identical - to the anomalous Hall 
conductivity. 

We computed the orbital magnetization as a function 
of the filling n by Eq. Q in the metallic region at t' = 0.5 
and in the insulating region at t' = 1.0. The behaviour of 
the orbital magnetization as the Fermi energy varies from 
the bottom of the lowest band to the top of the second 
band is displayed in Fig. [3J The orbital magnetization 
in the insulating phase at t' = 1 shows a discontinuity 
at n = 0.25, because the integral of the Pancharatnam- 
Berry curvature over the Brillouin zone is non-zero and 
quantized. The anomalous quantized Hall conductivity is 
the ratio of the discontinuity in the orbital magnetization 
to Eg I (2e), where E g is the gap between the first and the 
second band. 
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FIG. 3. (Color online) Orbital magnetization of the model as 
a function of the filling n for £'=0.5,1.0 

Thus far, the calculations we have presented were ob- 
tained in the non- interacting limit (U = 0). We now 
study the model as a function of the Hubbard interac- 
tion U to determine the region in the parameter space 
(t f — U) where the topological effects persist. The defini- 
tion of the Chern invariant can be generalized to interact- 
ing systems in terms of the one-particle Green's function 
(see, e.g., Eq. (66) of Ref. [18). It is expected to remain 
invariant provided the interaction does not close the gap 
and introduce new poles in the Green's function. It will 
also manifest itself as chiral edge states in the interacting 
system. We therefore address two questions: What is the 
region of parameter space where the gap persists? Are 
there chiral edge states in that region? 



We study the interacting theory using Cluster Pertur- 
bation Theory [11] (CPT) and the Variational Cluster 
Approximation [12]. See Supplemental Material for a 
brief review of the technique. 
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FIG. 4. (Color online) Clusters used in this work. Left: 6-site 
cluster used for the two-dimensional model. Right: 10-site 
cluster used for one-dimensional zigzag ribbons. Intercluster 
links are represented by dashed lines and the repeated unit is 
shaded in gray. 

We have applied the VCA to the current system 
by treating the cluster spin magnetization M c and the 
cluster chemical potential \i c as variational parameters. 
Fig. [4] illustrates the clusters that were used in this work. 
Allowing ji c to be different from ji and adopting the value 
that makes the Potthoff functional ^(M c ,/i c ) stationary 
ensures thermodynamic consistency, i.e., that the elec- 
tron density n calculated from the CPT Green function 
coincides with —dVt/dn [19]. This procedure allows us to 
calculate a more accurate value of the electron density n, 
compared with the simple CPT result. 




FIG. 5. (Color online). The phase diagram of the model in 
U-t' plane. 

We have scanned several values of the interaction U 
and of the spin-orbit coupling t' : in order to find whether 
the system remains gapped at quarter- filling. Fig. [5] 
shows the phase diagram thus obtained, on the t f — U 
plane. The curve shown is the phase boundary, i.e., the 
critical value t' c at which the Fermi energy moves from the 
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gap to within the band, as a function of U. In the insulat- 
ing region (right) the system remains in the quantum Hall 
state. Left of the line, the gap disappears and the sys- 
tem becomes metallic. We call this a chiral metal because 
non quantized Hall current flows along its boundary. We 
conclude that the topology of the band is protected for 
U > 0, and that a gapped state exists for t' > 0.5 if U is 
large enough. 
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FIG. 6. Top: Single particle spectral function, as a function of 
uj, for spin-up electrons and wavevectors along high-symmetry 
directions at t' — 1 and U — 0,4. Bottom: the same, at 
t' = 2/3. In all cases the system is quarter filled. 
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FIG. 7. (Color online) Top: Single particle spectral function 
of zigzag ribbons of 10 sites, as a function of energy uj. Top: 
up spins; bottome: down spins. 

tional anomalous quantum Hall states in this model when 
the bands are partially filled. Work in this direction is in 
progress. 

We are grateful to G. Baskaran, H.R. Krishnamurthy 
and A.-M.S. Tremblay for useful discussions. Computa- 
tional resources were provided by Compute Canada and 
Calcul Quebec. 



Fig. [6] displays the single-particle spectral weight 
A^(k,uj) of up-spins, as a function of frequency for 
wavevectors k along high-symmetry directions. Three 
of the four plots are taken from the gapped phase, and 
one (bottom left) from the metallic phase. Although the 
gap is small at U = 4 and t' = 2/3 (bottom right), it 
is still nonzero. A Lorenzian broadening of r] = 0.03t 
has been used to make these plots, but the phase dia- 
gram itself was obtained by using several values of r] and 
extrapolating the density of states towards r] = 0. 

We also calculated the spectral weight of zizgag rib- 
bons of width 10 at t' = 1 and quarter filling. The system 
is shown on Fig. QB) and the spectral weight on Fig. [7| 
The bottom four bands are nearly degenerate at k = tt 
in the presence of interactions and are still clearly sepa- 
rated from the chiral edge state (shown by the red arrow) 
near that wavevector. This indicates that the chiral edge 
state persists at U > 0: the interaction does not destroy 
the topoplogy of the bands. 

In conclusion, we have shown that time reversal break- 
ing spin-orbit interactions, in a model that can be real- 
ized in fermionic cold atom systems [2H4] , lead to anoma- 
lous quantum Hall states at - and | filling. These states 
are robust against interactions and have a clear experi- 
mental signature, if the orbital angular momentum can 
be measured. Our result motivates the search for frac- 
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SUPPLEMENTARY MATERIAL 



The Chern numbers 



Diagonalisation at U = 

In this section we describe the details of the diagonal- 
ization of the model in the U = limit. We label sites 
by (21,^2 5^0 where a = 1,2 labels the sublattice. The 
unit cell positions span a triangular lattice and are given 
by ziei + ^2^2, where the basis vectors are taken to be 

ei(2) = (T§,^). 



We define the Fourier tranform 



; ,k-(iie 1 +i 2 e 2 ) f 



(5) 



Let P a = |(1 + tV a ) and fci(2) = k * ^1(2) • The Hamilto- 
nian in momentum space can then be written as {t = 1), 



^-E( c li c k >2 ) f E t (k ) 



E(k) \ / ck,i 



M c k , 



(6) 



where E(k) = P 3 + Pie ifel + P 2 e~ ik2 . We define a x and 
a y to be the two Pauli matrices with sublattice indices. 
Further defining a* = (a x ± ia v )/2 the single particle 
hamiltonian can be written as, 



h(k) = a+® E(k) + a~ ® S f (k) 



(7) 



The single-particle spectrum is completely determined in 
terms of the spectrum of the positive semi-definite matrix 



^ = f*f+^(t') 2 +'-B-a 



2" 



(8) 



where 



/=2(l + e' 



2 ) 



and 



B = (1 — f' sin ki + cos &2 + cos £3) bi 
+ (1 + cos ki — t' sin fc 2 + cos ^3) t>2 
+ (1 + cos k\ + cos &2 — sin ^3) b3 (9) 

If ^> ± (k) are the eigenvectors of B(k) -cr, with eigenvalues 
±|B(k)|, then they are also the eigenvectors of E(k)Et(k) 
with eigenvalues, 



t' 



4(k) = r/+ 7 (0 2 ± oiB(k)i 



The four-component vectors 



W " 72 I ±^(-k) 



(10) 



(11) 



are then the eigenvectors of the single-particle Hamilto- 
nian defined by Eq. (pi) with eigenvalues ±e±(k). 



The Pancharatnam-Berry curvature for each band is 
given by 

PB pp '(k) = |J ^ pp/ (k^dj^ pp \k) - H.c.) (12) 



From Eq. (11) it is clear that 



PB^(k)=p^(Kk) + K-k)) 



6(k) 



8tt 



B(k) -^B(k) x ^-B(k) (13) 



It is not difficult to prove that |B| 7^ for all points in the 
Brillouin zone, including the Dirac points, for any non- 
zero value of t' . Thus PB PP (k) is well defined through- 
out the Brillouin zone. We have numerically computed 
J k 6(k) to be 1. Thus the Chern number of band pp' is p ': 
The highest and lowest energy bands, e±+, have Chern 
numbers +1 whereas the middle two bands, e±_, have 
Chern numbers —1. 

In the presence of interactions, the Chern number may 
still be calculated in principle, using the following expres- 
sion [18] : 

N= —^2 J d 3 ke^ x tv {Gd^G^GdM^GdxG- 1 } 

(14) 

where G stands for the Green function, the integral is 
taken over frequency and momentum and the trace is 
taken over spin and band indices. Greek indices are 
space-time indices running from to 2. In the non- 
interacting case, this expression coincides with the inte- 
gral of the Pancharatnam-Berry curvature over occupied 
states of the Brillouin zone. In the interacting case, the 
Green function will be deformed, but as long as a gap 
persists at the Fermi level, the Chern number should not 
be affected, as it is robust against smooth deformations 
that do not introduce or remove any poles in G. 



CPT and VCA 

CPT is an approximation scheme for the one-electron 
Green function G(gj) within Hubbard-like models. [TTJ[20j 
|2T] It proceeds by dividing the infinite lattice 7 into a su- 
perlattice T of identical clusters of L sites each. The lat- 
tice Hamiltonian H is written as H = H c +Ht, where H c 
is the cluster Hamiltonian, obtained by severing the hop- 
ping terms between different clusters, which are put into 
Ht- If T is the matrix of inter-cluster hopping terms and 
G c (uj) the exact Green function of the cluster. Because 
of the periodicity of the superlattice, T can be expressed 
as a function of the reduced wavevector k and as a matrix 
in site indices within the cluster: T a b(k). Likewise, G c is 
a matrix in cluster site indices only, since all clusters are 
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identical: G c ah {uo). Thus, hopping matrices and Green 
functions in what follows will be k-dependent matrices 
of order L, the number of sites within each cluster. The 
fundamental result of CPT is 

G-^cj) = G c -\uj) - T(k) (15) 



dft/dfi = — n gives us a reliable estimate of the electron 
density. VCA provides estimates of order parameters, 
much like mean-field theory, but is quite superior to it 
because the Hamiltonian remains fully interacting (no 
factorization of the interaction) and spatial correlations 
are treated exactly within the cluster. 



In practice G c (u) is calculated numerically by the Lanc- 
zos method and the cluster must be small enough for this 
to be possible. Because the lattice tiling breaks the orig- 
inal translation invariance of the lattice, a prescription 
is needed to restore the translation invariance of the re- 
sulting Green function. The CPT prescription for this 
periodization is 



a, 6 



(16) 



where now k belongs to the Brillouin zone of the original 
lattice. This formula is exact in both the strong (t — >• 0) 
and the weak (U —> 0) coupling limits. 

Once the approximate interacting Green function can 
be calculated, various quantities can be calculated, such 
as the electron density n(/i) as a function of chemical 
potential, or the spectral function A(k, uo) and its integral 
over wave vectors, the density of states N(uS). 

The Variational Cluster Approximation (VCA) is an 
extension of CPT in which parameters of the cluster 
Hamiltonian H c may be treated variationally, accord- 
ing to Potthoff's Self-Energy Functional Theory (SFT) 
[I2j[22]. In particular, it allows the emergence of sponta- 
neously broken symmetries and provides an approximate 
value for the system's grand potential Q. Technically, 
VCA proceeds by minimizing the following quantity: 



fi(ft) = fi c (ft)- 



V(2^? lnd6t i ~ T ( k ) G ( k '* w ) 



(17) 

where fl c (h) is the grand potential of the cluster alone 
(obtained in the exact diagonalization process). The in- 
tegral over frequencies is carried over the positive imag- 
inary axis, and h denotes collectively the parameters of 
the cluster Hamiltonian H c that are treated variation- 
ally; these must be the coefficients of one-body operators. 
At the optimal value h* , Q(h*) is the best estimate of 
the system's grand potential; in particular, its derivative 



[3; 



[5] 

[6] 
[7] 
[8] 
[9] 



no 



[ii 

[12; 
[13; 



[14 



A. Kitaev, Annals of Physics, 321, 2 (2006) 



L.-M. Duan, E . Demler, and M. D. Lukin, Phys. Rev. 
Lett., 91, 090402 (2003)| 



C. Zhang, V. W. Scarola, S. Tewari, and S. Das Sarma, 
Proceedings of the National Academy of Sciences, 104, 
18415 (2007)T 

R. Jordens, N. Strohmaier, K. Giinter, H. Moritz, and 
T. Esslinger, Nature, 455, 204 (2008). 

D. J. Thoul ess, M. Kohmoto, M. P . Nightingale, and 
M. den Nijs JPhys. Rev. Lett., 49, 405 (1982)] 

Y. Hatsugai, lPhys. Rev. Lett., 71, 3697 (1993)| 

F. D. M. Haldane Phys. Rev. Lett., 61, 2015 (1988). 

F. D. M. Haldane, Phy s. Rev Lett., 9 3, 20660 2 (2004)| 

M. Z. H asan and C. L. Kane, |Rev. Mod. Phys., 82, 3045] 

(2010)| 

X.-L. Q i and S.-C. Zhang, Rev. Mod. Phys., 83, 1057 



(2011 

D. Senechal, D. Perez, and M. Pioro-Ladriere, Phys. 
Rev. Lett., 84, 522 (2000). 
M. Potthoff, Eur. Phys. J. B, 32, 429 (2003). 
J. Shi, G. Vignale, D. X iao, and Q. Niu, |Phys. Rev" 
Lett., 99, 197202 (2007")] 



D. Xiao, M.-C. Chang, and Q. Niu, Rev. Mod. Phys., 



82, 1959 (2010) 

D. Xiao , J. Shi, and Q. Niu, Phys. Rev. Lett., 95, 137204 
(2005) 

T. Th onhauser, D. Cere soli, D. V anderbilt, and R. Resta, 
|Phys. Rev. Lett., 95, 137205 (2005)] 



D. Xiao, Y. Yao, Z. Fang, and Q. Niu, Phys. Rev. Lett., 
97, 026603 (2006)| 

G. Volovik, Physics Reports, 351, 195 (2001). 
M. Aichhorn, E. Arrigoni, M. Potthoff, and W. Hanke, 



[15 

w 

lir 

[18 
[19 

|Phys. Rev. B, 74, 235117 (2006")] 
[20] D. Senechal, D. Perez, and D. Plouffe, Phys. Rev. B, 66, 
075129 (2002). 

[21] D. Senechal, in Theoretical methods for Strongly Cor- 
related Systems, Springer Series in Solid-State Sciences, 
Vol. 171, edited by A. Avella and F. Mancini (Springer, 
2012) Chap. 8, pp. 237-269. 

[22] M. Potthoff, in Theoretical methods for Strongly Cor- 
related Systems, Springer Series in Solid-State Sciences, 
Vol. 171, edited by A. Avella and F. Mancini (Springer, 
2012) Chap. 9. 



